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Abstract 

We present an unrestricted Hartree-Fock computation of charge-ordering instabilities of two-dimensional 
metals with antiferromagnetic exchange interactions, allowing for arbitrary ordering wavevectors and in- 
ternal wavefunctions of the particle-hole pair condensate. We find that the ordering has a dominant d 
symmetry of rotations about lattice points for a range of ordering wavevectors, including those observed 
in recent experiments at low temperatures on YBa2Cu30 y . This d symmetry implies the charge ordering 
is primarily on the bonds of the Cu lattice, and we propose incommensurate bond order parameters for the 
underdoped cuprates. The field theory for the onset of Neel order in a metal has an emergent pseudospin 
symmetry which 'rotates' <i-wave Cooper pairs to particle-hole pairs (Metlitski et ah, Phys. Rev. B 82, 
075128 (2010)): our results show that this symmetry has consequences even when the spin correlations are 
short-ranged and incommensurate. 
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A remarkable series of experiments [fTl-[T3ll have shed new light on the underdoped region of 
the cuprate high temperature superconductors. These experiments detect a bi-directional charge 
density wave with a period of about 4 lattice spacings at low hole densities and low temperatures 
(T). This order is co-incident with regions of the phase diagram where quantum oscillations 
[fT4| were observed in YBa2Cu 3 3 „ strongly supporting the hypothesis 1171 [T514T81I that the charge 
density wave is responsible for the Fermi pockets leading to quantum oscillations. Some of the 
experiments [|3l |5]-[H [EH EE2l [191 I2Q1 also indicate that there is negligible modulation of the charge 
density on the Cu sites; instead, it is primarily a bond density wave, with modulations in spin- 
singlet observables on the Cu-Cu links, such as the electron kinetic energy. 

This paper presents a Hartree-Fock computation of charge-ordering instabilities of a two- 
dimensional metal of electrons with antiferromagnetic exchange interactions (described by a l t-F 
model). We allow the charge-ordering to appear at any wavevector, Q, and also allow an arbitrary 
internal wavefunction, Ag(fc) for the spin-singlet particle-hole pair condensate responsible for the 
density wave order (here Q is the center-of-mass momentum of the particle-hole pair, and k is the 
relative momentum). We show that this freedom leads to significant insight, despite the simplicity 
of our method. We find that for a range of small Q (more precisely, in the "7~ preserved' region 
of Fig. [2]), including those observed so far in the experiments ffTT-461 [814T2TI at low T, the dominant 
structure of the internal wavefunction has a rf-wave form lETI . with Aq(&) ~ (cos k x - cos k y ) for 
a band- structure appropriate for the cuprates. This d symmetry implies that the charge order is 
located primarily on the Cu-Cu links, there is little modulation of the charge density on the Cu 
sites, and time-reversal symmetry (7~) is preserved. We refer to this class of charge order as an 
'incommensurate J- wave bond order'. Our computation also allows for other spin-singlet orders, 
such as Ising-nematic [|22l - [24l . 'd-density wave' ll25l - [27ll . and 'circulating currents' 11281 . the last 
two of which break 7~: they are all less preferred than the incommensurate d- wave bond order in 
the underdoped region, while Ising-nematic order is preferred at larger doping 11291 . 

The preferred value of Q in our Hartree-Fock computation in a metal has the form Q = 
(+Q,„, +Qm) [21J; similar orders have appeared in recent computations |[2"Tll2"9"ll30l using the renor- 
malization group and other methods. At low doping, we find that Q m « Q , where Qo is defined 
geometrically from the 'hot spots' on the Fermi surface, as shown in Fig. [T] (see also Fig. [3] for a 
comparison between the values of Q m and Q ). Remarkably, the hot spots of commensurate Neel 
order play a crucial role when the antiferromagnetic correlations are short-ranged, and even when 
they are incommensurate. Recent field-theoretic studies ||2T1 [301 focused on the Fermi surface in 
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FIG. 1 : Fermi surface of the hole-doped cuprates, showing the Brillouin zone boundary to antiferromag- 
netism at K = (n,7t) (dashed lines), the hot spots (filled circles), and the wavevectors (Qo,0) and (<2o> £2o)- 

the immediate vicinity of these hot spots, and this connection allows us to interpret our Hartree- 
Fock results in terms of an emergent approximate pseudospin symmetry I12H . The pseudospin 
symmetry rotates 'rotates' J- wave Cooper pairs to particle-hole pairs: the Cooper pair amplitude 
A s (k) rotates into Ag(fc), which explains the predominant d symmetry of the latter. Our results 
show that the pseudospin symmetry is a good guide to picking optimum states in lattice computa- 
tions on models with short spin correlation lengths, even though the symmetry is exact only in a 
continuum limit where the spin correlation length becomes very large. 

As in Ref. [|2Tl[30ll . we propose that the high T pseudogap is a metal with a fluctuating multi- 
dimensional order with both a superconducting component, A s (k), and a bond order component 
Aq(fc) over a range of small Q. At lower T, superconductivity appears by the polarization of 
this order along A s (k). Subsequent static charge-ordering requires computation of the non-zero 
Q instabilities within the superconductor. Fortunately, the latter computation has already been 
performed in closely related models PT1I321I : bond order modulations were found with Q along 
the (1,0), (0, 1) directions, as is the case in observations at low T [Q3-|6l l8Ul2ll. In our metallic 
computations, there is a 'valley' of stability from (Q m , Q m ) to (Q m , 0) and (0, Q m ), but the global 
minimum is at (Q m , Q m ) (see Fig. [2]); within the superconducting phase in zero magnetic field, the 
balance can evidently be tipped in favor of bond order near (+Q m ,0) and (0, ±Q m ). The choice 
of ordering wavevectors between (+Q m , +Q,„) and (+Q m , 0), (0, +Q,„) surely depends upon details 
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of the Hamiltonian, the value of T, and the presence of a magnetic field, and is perhaps not ac- 
curately estimated by our present simple Hartree-Fock computation. Nevertheless, we expect the 
predominant d symmetry of particle -hole pair condensate, Ag(fc), to be robust for \Q\ < 2Q m , for 
the same reason it is robust for the particle-particle condensate, A s (k), of the superconductor. 

We examine the following Hamiltonian of electrons on a square lattice of sites i at positions r, 
with annihilation operators c ia , where a =T, I is a spin label: 



(1) 



Here n is the chemical potential, t t j are the electron hopping amplitudes, Ju are exchange interac- 
tions, and the electron spin operator §i = |ct,(? Qs gC^ 8 , with a the Pauli matrices. The pseudospin 
symmetry acts as on the Nambu spinor *P !Q , = (c ia , e a pc]^) as a SU(2) rotation V, in Nambu space 
under which — > V,^. A key property is that 5*, = l^d^Y^, and this is invariant under 
the pseudospin transformation. Consequently the exchange interaction is invariant under indepen- 
dent rotations V, on each lattice site [33 J, and this gauge invariance was exploited in the study of 
spin liquid ground states of Mott insulators Il27l . The pseudospin symmetry is completely broken 
by the f, 7 terms in H, and so it was expected that pseudospin symmetry plays no role in metallic 
states, except those that are proximate to certain spin liquids [IT71 |3"4l . Here, we are interested 
in metallic states proximate to systems with long-range antiferromagnetism; surprisingly, it was 
shown in Ref. ETTl . that an analog of the pseudospin gauge symmetry of Refs. 11331 reappears in 
the critical theory of the antiferromagnetic quantum critical point in a conventional metal, as 4 
independent global SU(2) pseudospin rotations, one for each pair of hot spots. These rotations 
serve to map the d- wave Cooper pairing A s (k) to the rf-wave bond order Aq(fc), as is also evident 
in our computations below. 

Note that H does not contain an explicit on-site interaction, the 'Hubbard U' . Both the Cooper 
pair and the bond order have small on-site components because of the d symmetry, and so U is not 
important in selecting the ordering instabilities. The effects of U can be accounted for by 'slave 
particle' methods ll25l I3TI . and its main consequence is a renormalization of the quasiparticle 
dispersion. Finally, such local interactions are irrelevant in the field theory of Ref. I12TTI . 

For our charge-ordering Hartree-Fock analysis, we need the best variational estimate for the 
mean-field Hamiltonian 

H M f = J] (-friiJ ~ f iJ ~ A o) c l c ja (2) 
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where the non-local charge order A, ; is written as 



V 

Q L fc 



e^W)/ 2 , (3) 



with V the system volume. This expression highlights the physical interpretation of Ag(fc): (z) if 
Aq(fc) is a constant independent of k (i.e. s-wave) then we have an ordinary site charge density 
wave at wavevector Q with only A,,- non-zero; (ii) if Aq(fc) ~ c\ cos k x + c 2 cos k y (d- and extended 
s-wave) then we have bond order at wavevector Q with A, 7 non-zero only if i and j are nearest 
neighbors. Also note that hermiticity requires Ag(fc) = A_q(A;), and under time-reversal T : 
A Q (fc) -> A Q (-fc). 

All the functions Aq(fc) are variational parameters, to be optimized by minimizing the free 
energy by F < F MF + (H - H MF ) MF , where the average is over a thermal ensemble defined by 
H M f- Here, we expand the r.h.s. in powers of Aq(A;), and replace the inequality by an equality. 
To quadratic order in Aq, we write the result in terms of hermitian functional operators on the 
Brillouin zone as 

f= Yj A b( fc ) A /n Q (fc)M Q (fc, fe') ^n Q {k')k Q (k') + ... (4) 

k,k',Q 

where the kernel is 



M Q (k, fe') = 6 kM + - X o(k - k') ^n Q (fc)n Q (fc') (5) 

while the polarizability and susceptibility are 

rr ru, mk + Q/2))-f(s(k-Q/2)) 1 v iq , } 



(6) 



j 

with s(k) the electron dispersion associated with fy, and / the Fermi function. From Eq. (|4]) we see 
that the linear charge-ordering instability of the metal occurs via condensation in the eigenmodes 
of the operator Mq(k, k') with the lowest eigenvalues. We have chosen the specific forms of the 
kernel in Eq. ([5]) so that we need only solve the following eigenvalue problem 

£ Jn Q (k)xo(k - k') ^(fcO0 Q (fc') = A Q( p Q {k) 

k' 

for the minimum eigenvalues Aq and corresponding eigenvectors (f>Q(k), and their structure is 
independent of the overall strength of the interaction xo- The charge-order will then be Aq(A0 oc 
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4>Q(k)/ ^TlQ(k). Our principal numerical results below are on the Q dependence of Aq, and on 
the k dependence of Ag(fc) so obtained. 

We also solved for the corresponding instability of the metal to the superconductor. In this 
case H M f has the charge-ordering term A, 7 replaced by the pairing term - Yik As(fc)cfcT c -fcl + H.c, 
and the subsequent expressions have the replacements Aq(&) — » A s (k), Mct(k, k') — » Ms(k, k'), 
n e (fc) —> Tls(k), Aq —¥ A s , with H s (k) = (1 - 2/(e(fe)))/(2e(fc)). In particular, the expression for 
the kernel Ms(k, k') in terms of II s (fc) has a form identical to Eq. (|5]), a key consequence of the 
pseudospin symmetry of the exchange interaction. Note also that for dispersions with e(k + Q) = 
-e(k) we have IIq = H s and so Mq = Ms ; Ref. [|2T| pointed out that the dispersion obeys such 
a relationship holds close to the hot spots of a generic Fermi surface for Q = (+Qo, ±Qo) (see 
Fig.[T])), and this then establishes the pseudospin rotation symmetry between A s and Aq(&)- 

We assume an electronic dispersion s(k) = -2t\(cos(k x ) + cos(k } ^J - 4t 2 cos(k x )cos(k y ) - 
2t 3 (cos(2k x ) + cos(2^,)) ~M an d a susceptibility Xo(q) which is peaked near the antiferromagnetic 
wavevector 

- V A 

X ° {q) ~ V 4(r 2 + 2(2 - cos(q x - K x ) - cos(q y - K y ))) ' (7) 

where ^ is the antiferromagnetic correlation length, the sum extends over K = n{\ - 5)), + 
(n(l-6), n), and we used both the commensurate case 8 = and the incommensurate case 8 = 1/4, 
with little difference between the results. We only need a short spin correlation length, and 
indeed obtained very similar results even for the case where Xo(q) was obtained from Eq. ([6]) with 
only a nearest-neighbor Jy. We diagonalized the kernels after discretizing the Brillouin zone to L 2 
points with L up to 80, and the results below are for t\ = l,t 2 = -0.32, and f 3 = 0.128 for a range 
of values of T, fi, and 

Numerical results. For the full range of parameters examined, we consistently found that As 
was the minimal eigenvalue (indeed, BCS theory implies —As diverges logarithmically as T — » 0), 
and the corresponding eigenvector A s (fc) was well approximated by the rf-wave form ~ (cos k x - 
cos k y ) (see Table [Ij). So d-wave superconductivity is the primary instability. 

For the charge ordering instabilities, we show the Q dependence of Aq in Figs [2] and in the 
supplement. We characterize the corresponding eigenvectors A S g(fc) using orthonormal basis 
functions, iff y (k) of the square lattice space group: 

A Q {k) = Y J CQ,y^y{k) (8) 

r 

where cq^ are numerical coefficients collected in Table |IJ Depending upon the symmetry of Q 
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T broken 




FIG. 2: Plot of Aq/A, where Aq is the smallest charge order eigenvalue, as a function of Q x and Q y . We 
used ju = -1.11856, £ = 2, T = 0.06, 6=1/4 and L = 64. Charge order appears when Aq < -1, which 
happens when A is large enough. The global minimum is at (Q m , Q m ) and Q m is plotted in Fig. [3] as a 
function of /u. Notice also the blue valleys extending from (Q m , Q m ) to (Q m , 0) and (0, Q m ). The region with 
time-reversal, T, preserved has the eigenfunctions Ag(-fc) = Ag(fc) which are predominantly d, while the 
region with T broken has Aq(-A;) = -Ag(fc), as shown for some values of Q in TablejlJ 

(in particular, the little group of the wavevector Q) and of the eigenvector, some of the cq 7 may 
be exactly zero. But for a generic Q, only time-reversal constrains the values of cq 7 , and we are 
allowed to have an admixture of many basis functions. Nevertheless, only a small number of basis 
functions have appreciable coefficients, and so Eq. ^ represents a useful expansion. 

The global minimum of Aq is at a wavevector along the diagonal with Q = (Q m , Q m ), and we 
show a plot of Q m as a function of chemical potential in Fig. [3] We also show the corresponding 
values of Q as defined in Fig. [TJ for small doping we see that Q m « Q , one of our key results: 
the pseudospin symmetry of the hot-spot theory of Ref. ETTl is a good guide to determining the 
ordering even for models with short-range, incommensurate, antiferromagnetic spin correlations. 
At larger doping, after the chemical potential crosses the van-Hove singularity il29ll . there are no 
hot spots, and we find Q m = 0. For Q = (Q m , Q m ), Table [j] shows that that Ag(fc) is predominantly 
d, with a small admixture of g. For Q = (Q m ,0), Aq remains predominantly d, but now has a 
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tffy(k) 


Q = 


Q = 


Q = 


Q = 


As(fc) 






(.Qtni Qm) 


(Qm,0) 


(0,0) 


(n,n) 
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1 





-0.226 











s' 


cos k x + cos jL 

y 





0.040 











s" 


cos(2^ x ) + cos(2£y) 





-0.051 
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cos k x - cos k v 


0.993 


0.964 


0.997 





0.998 


d' 


cos(2^) - cos(2fc y ) 


- 0.058 


-0.057 


-0.044 





-0.047 
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2 sin/c,- sin k v 

- x y 

















Px 


V2 sin k x 











0.706 





Py 


V2 sin & y 











-0.706 





g 


(cos k x - cos fcy) 


-0.010 
















x V8 sin k x sin £ y 













TABLE I: Values of cq iT in the expansion for Ag(fc) in Eq. ([8]) for various values Q and y. The values 
of cq jT are normalized so that |cq, 7 | 2 - 1, where the sum over y includes the small contributions from 
higher order basis functions not shown above. Values shown as are constrained to be exactly zero by 
symmetry. The last column shows the coefficients in the corresponding expansion for A$(k). Parameters 
are as in Fig. [2} and Q m = AnjW. 

small s components (32). 

At Q = 0, we find that Ag(fc) is purely d: this corresponds to Ising-nematic order Il22] - l24l . The 
7~-breaking 'circulating-current' order of Ref. Il28l has a p x<y form for Aq(&) at Q = 0, but this 
does not appear as a lowest eigenvalue, and so is not present in Fig. [2] Finally, Aq also has a broad 
local minimum at Q = (n,n): here Ag(fc) does have the p X:) , form which breaks T, and leads to 
the state with spontaneous orbital currents [|25l - l27ll . 

Experiments IfTMol [8UT21 have observed charge ordering at Q = (Q m ,0), (0, Q m ) at low T. 
Choosing the largest 2 components at this wavevector from Table [IJ we have 

A s + A (/ (cos k x - cos ky) , Q = (+Q m , 0) 
A. s - A d (cos k x - cos k y ) , Q = (0, +Q m ) 
with Av/Aj = -0.234. Similarly, we can have bond-ordering along Q = (+Q m , +Q m ) with only 
Ad non-zero. We present implications of these orders for X-ray scattering, photoemission and 
scanning tunneling microscopy in the supplement. 



A Q (fc) 



(9) 
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FIG. 3: Plot of Q m (circles), where the minimum of Aq occurs at Q - (±Q m , ±Q m ). Also shown are the 
corresponding values of Qq (squares), as defined by the hotspots on the Fermi surface in Fig. [T] The near 
equality of Q m and <2o is evidence for the pseudospin symmetry; note that this holds even though Xo(o) 
in Eq. ^ is peaked at the wavevectors K = (n, ±37r/4), (±37r/4, n), as is the case in many hole-doped 
cuprates. 

Our evidence for pseudospin symmetry between Cooper pairing and charge order should have 
significant implications for the dynamics of these orders, which have been studied recently in 
Ref. 11351 . For the phase diagram of the hole-doped cuprates, our model has a T = quantum- 
critical point near optimal doping associated with disappearance of this bond order [JVl [3TJ . An 
important challenge is to use such a critical point to describe the evolution of the Fermi surface 
ifTTll . and the 'strange' metal. 
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SUPPLEMENTARY MATERIAL 



First, we give further details on the function Aq in Fig. [2] In Fig.|4j we plot Aq along different 
lines in the Brillouin zone, and also indicate the regions where T is preserved and broken. 

X Q /A 
-o.i 

-0.15 
-0.2 
-0.25 

0.25 0.5 0.75 ^ , 1 

FIG. 4: Plot of the eigenvalue of Fig. [2] along the Brillouin zone diagonal with Q = (Q, Q) (full line), along 
the line Q = (Q,„, Q) (dotted blue line), now with L = 80. The eigenfunction Aq(&) has predominant d 
symmetry (as in the state of Ref . EH ) with T preserved to the left of the filled circles, and predominant p x ^ y 
symmetry with 7~ broken (as in the state of Refs. E5ll26l ) to the right of the filled circles. The Q = (0, 0) 
point corresponds to Ising nematic order IF22TI241 . 

Next, we describe properties of the bond-ordered state in Eq. ([£]). Inserting Eq. ([£]) into Eq. ([3]), 
we see that the real space order parameter A !7 is non-zero only when i = j, or when i and j 
are nearest neighbors. The values of A it correspond to an ordinary on-site charge density wave 
on the Cu sites at wavevectors Q = (0, ±Q m ), (+Q m ,0) with amplitude proportional to A 4 . The 
larger component of the ordering is however the bond-density wave given by A, ; - with i,j nearest 
neighbors, whose amplitude is proportional to A^. We show plots of the values of A (J on the 
bonds of the square lattice in Fig. [5] and [6] Fig [5] contains the case of uni-directional order only 
at the wavevectors Q = (+Q m ,0), while Fig. [6] is the case of bi-directional order at wavevectors 
Q = (±Q m ,0) and Q = (0, ±Q m ). 

For completeness, we also show the corresponding plots for ordering along Q = (±Q m , ±Q m ) 
in Figs[7]and[8j these appeared earlier in Figs. 22 and 23 in Ref. j|2T| at a different period. Note that 
the difference between bi-directional order at Q = (+Q m ,0) and Q = (0, ±Q m ) in Fig. [6] and bi- 
directional order at Q = +(Q m , Q m ) and Q = +{Q m , -Q m ) in Fig.[8]is subtle, and not immediately 
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FIG. 5: The bond-component of the ordering in Eq. (|9]>, proportional to Ad, for the case of uni-directional 
order at Q = (±Q,„, 0) where Q m = AnjW. There is an additional site-component, proportional to A s , which 
is not shown. The lines intersect at the Cu sites, and the colored squares are on the O sites: the colors are 
therefore a measure of the charge density (or other spectral properties) on the O sites. 
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FIG. 6: As in Fig. [5] but for the case of bi-directional order at Q = (±Q m , 0) and Q = (0, ±Q m ). 
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FIG. 7: As in Fig.[5j but for the case of uni-directional order at Q = ±(Q m , Qm)- We have chosen Ag(fc) to 
be purely d, which is an excellent approximation to the state in Table [I] In this case Ay is non-zero only if i 
and j are nearest neighbors, and these are shown above; there is no density wave on the Cu sites. This plot 
also appeared in Ref. 11211 with a different period. 
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FIG. 8: As in Fig. [7] but for the case of bi-directional order at Q = ±(Q m , Q,„) and Q = ±(Q,„, -Qm)- 
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apparent at first glance: the periods along the x and y axes appear the same. However, the Fourier 
transforms of these two cases are distinct. 

The four plots in Fig. [5]j8] together contain information that should be useful in interpreting 
scanning tunneling microscopy and X-ray scattering experiments: the colors can be viewed as a 
measure of any observable on the O site which is invariant under time-reversal and spin rotation. 
Most simply, such an observable is the charge density on the O site, but any spectral property 
of the O atom also qualifies, and the latter can have readily measurable consequences in such 
experiments. 

Finally, we consider the electronic spectral function in the presence of bond-ordering in a metal. 
This is obtained by diagonalizing the following Hamiltonian 



H b = Yj <k)cl a c ka - J> Q (fc + Q/2)cl Q a c ka 

k L Q 



(10) 



where the sum over k extends over the complete Brillouin zone of the square lattice. For the case 
of bi-directional order in Eq. ([9]), the sum over Q extends over the 4 values (+Q m , 0) and (0, +Q m ). 



Some care must be taken in evaluating the wavevector Q/2 in the argument of Aq in Eq. (10) 
as it is not invariant under translation of Q by a reciprocal lattice vector of the square lattice: in 
each term, we take the momenta k and k + Q to be separated by exactly Q (and not modulo a 
reciprocal lattice vector), and then Aq(& + Q/2) is evaluated at the midpoint between them. For 
Q m = 4n/l 1, determining the spectrum of H h involves diagonalizing al21xl21 matrix for each k. 
From the eigenfunctions and eigenvectors we computed the imaginary part of the single-electron 
Green's function, ImG^ki^ + Z? 7X the quantity related to the photoemission spectrum. For the 
bi-directional ordering along Q = (+Q m , 0), (0, +Q m ) of Eq. ([9]) the result is shown in Fig. |9j The 
corresponding result for bi-directional ordering along Q = +{Q m , Q m ), +{Q m , -Qm) is in Fig. [TOj 
in this case A. v = by symmetry, and only Aj was non-zero. 

The stability of the Fermi arc in the 'nodal' region (k x « k y ) is enhanced 021 l36l because of 
the weak coupling to the charge order, arising from the predominant d character of Eq. ([9]). In the 
anti-nodal region, the parent Fermi surface has been gapped out by the bond order, but 'shadows' 
are apparent at wavevectors shifted by the charge order. However, these Fermi surfaces should 
be easily broadened by impurity-induced phase-shifts in the charge ordering, while protecting the 
nodal arcs. Furthermore, contributions from the superconducting component of the pairing order 
parameter should also help fully gap out the antinodal region. 
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FIG. 9: Electron spectral density in the phase with bidirectional charge order at Q = (Q m ,0) and (0, Q,„) 
with Q m = An The left panel show \mGk,k{^> + iv) at <x» = and tj = 0.02; the right panel shows 
log [lmGk k(a> + irj)] for the same parameters, as a way of enhancing the low intensities. The dashed line 
is the underlying Fermi surface of the metal without charge order. The charge order is as in Eqs. (3]9) with 
Arf = 0.3, A s / Ad = -0.234, and other parameters as in Fig. [2] 




FIG. 10: As in Fig. [9], but for the case of bi-directional ordering along Q = ±(Q m , Q m ), ±(Q m , -Qm)- The 
charge order is as in Eqs. (|3|9[) with Aj = 0.3, and A s = is required by symmetry. 
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